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Abstract — Segmentation and volume measurement of the 
cardiac ventricles is an important issue in cardiac disease 
diagnosis and function assessment. Cardiac Magnetic 
Resonance Imaging (CMRI) is the current reference standard 
for the assessment of both left and right ventricular volumes 
and mass. Several methods have been proposed for 
segmentation and measurement of cardiac volumes like 
deformable models, active appearance, shape models, atlas 
based methods, etc. In this paper a novel method is proposed 
based on a parametric superellipse model fitting for 
segmentation and measurement volume of the left ventricle 
on cardiac Cine MRI images. Superellipses can be used to 
represent in a large variety of shapes. For fitting superellipse 
on MR images, a set of data points have been needed as a 
partial data. This data points resulted from a semi-automatic 
region growing method that segment the homogenous region 
of the left ventricle. Because of ellipsoid shape of left ventricle, 
fitting superellipse on cardiac cine MRI images has excellent 
accuracy. The results show better fitting and also less 
computation and time consuming compared to active contour 
methods, which is commonly used method for left ventricle 
segmentation. 

Index Terms — Cardiac MRI, Volume measurement, 
Deformable superellipse, Levenberg Marquardt algorithm, 
Region Growing, Gray space map 

I. Introduction 

Cardiovascular disease is a major cause of death in the 
western world, which is responsible for nearly half of all 
deaths in Europe and approximately a third of all deaths in 
USA. There are several modalities for diagnosing 
cardiovascular diseases such as magnetic resonance imaging 
(MRI), computed tomography (CT), angiography and etc. 
Many radiologists investigate MRI images, searching for 
representative projection of the heart and its chambers (right 
and left ventricles and atriums). 

The computation of clinical parameters to assess the 
cardiac function requires segmenting the cardiac ventricles. 
As the heart is a moving organ, images are acquired 
throughout the whole cardiac cycle, but two precise instants 
are of particular interest for the clinician: the time of maximum 
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filling, when the heart is the most dilated (end diastole, ED), 
and the time of greatest contraction (end systole, ES). To 
estimate these volumes it is necessary to apply invasive or 
non-invasive medical investigation techniques. The 
segmentation of ED and ES images of the Right Ventricles 
(RV) is currently performed manually in clinical routine. This 
long and tedious task, prone to intra- and inter-expert 
variability, requires about 20 min per ventricle by an expert 
clinician. The great need for automated methods has led to 
the development of a wide variety of segmentation methods 
[1]. Most of these methods compute a pixel wise 
correspondence between the current image (or frame) and 
model distributions of photometric (intensity based) and 
geometric properties of the target objects. In this area, 
methods can be categorized in thresholding [2], pixel 
classification [3-4], deformable models and atlas based 
segmentation. Among these methods deformable models have 
been greatly used as their flexibility, especially for this 
application [5-7], either on the form of 2D active contours or 
3D deformable surfaces. A review of papers on deformable 
models can be found in [8]. One of the disadvantages of 
these methods are being computationally expensive [9, 10] 
and reinitialization problem, that make them time consuming. 
Nonhomogenity region and missing boundary is other 
problems that affect result efficiency. Several works have 
been performed in literature in left ventricle volume 
measurement area. One of these methods is Active Shape 
Model (ASM) and Active Appearance Model (AAM). ASM 
consist of a statistical shape model, called Point Distribution 
Model (PDM), obtained by a PC A on the set of aligned shapes, 
and a method for searching the model in an image [11. 
Segmentation is performed by placing the model on the image, 
and iteratively estimating parameters using least square 
estimation. ASM have been extended to gray level modeling, 
yielding Active appearance models (AAM) [12], that represent 
both the shape and texture variability seen in a training set. 
This technique ensures to have a realistic solution, since 
only shapes similar to the training set are allowed. In [ 1 3 ] , the 
authors introduce a multistage hybrid model, arguing that 
AAM are optimized on global appearance but provide 
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imprecise border locations, whereas ASM have a great abil- 
ity to find local structures. They thus propose to concat- 
enate several independent matching phases, starting with an 
AAM early stage that positions the model onto the heart, 
followed by a hybrid ASM/AAM stage that allows for posi- 
tion refinement. A final stage of AAM aids in escaping a 
possible local minimum found during the ASM/AAM stage. 
Although this method provides accurate results, and, for the 
first time with AAM, results on the RV, the current model 
training has been limited to mid-ventricular, end diastolic im- 
ages. Estimation models have been done by Dulce et al [14]. 
They compared magnetic resonance data sets by using a 
modified Simpson-Rule on a biplane ellipsoid model and the 
3D data set of MR images. Goshtasby et al. [15] applied a 
two-stage algorithm for extraction of the ventricular cham- 
bers in flow-enhanced MR images. They approximate loca- 
tion and size of endocardialy surfaces by intensity 
thresholding and reposition points on the approximated sur- 
faces to nearest local gradient maxima. Then they fit a cylin- 
der into the point set. Weng et al. [16] proposed a learning 
based ventricle detection from MR and CT images based on 
a likelihood measure for region-of-interest detection. The 
contributions regarding ventricle segmentation using deform- 
able models have mainly dealt with the design of the data- 
driven term. Initially gradient based [17] and thus sensitive 
to noise, region-based terms have been introduced, that is 
based on a measure of region homogeneity [18]. In [19], the 
authors take into account the intensity distribution overlap 
that exists between myocardium and cavity, and background 
and myocardium, and introduce a new term in the functional 
that measures how close the overlaps are to a segmentation 
model, manually obtained in the first frame. A comprehensive 
review of segmentation methods in short axis cardiac MR 
images can be found in [20]. In this paper we used a paramet- 
ric superellipse model, for fitting on left ventricle. 
Superellipses are a flexible representation that naturally gen- 
eralizes ellipses. They can model a large variety of natural 
shapes, including ellipse, rectangles, parallelograms, and 
pinched diamonds, by changing a small number of param- 
eters [21]. Hence, the superellipse is a flexible primitive in 
machine vision and image processing. In addition, the 
squareness of the superellipse is similar to the eccentricity of 
an ellipse, which is very useful for pattern recognition. 
Superellipses were first formulated by Gardiner [22] . The three- 
dimensional superellipsoid version was popularized in com- 
puter graphics by Barr [23] and in computer vision by Pentland 
[24]. Several approaches have been suggested to determine 
the parameters of superellipses. In its original form, a 
superellipse is fitted using the least mean square (LMS) 
method, and the parameters are determined by minimizing 



the errors between the curve and the real data. Orthogonal 
and algebraic distances are used to calculate the errors [25]. 
When only partial data (data points) is available the param- 
eter estimation and Fitting superellipse to the data become 
unreliable [26]. There are several methods that can provide 
partial data for fitting superellipse on them. Some of these 
methods are manually and need user interaction. But the ad- 
vanced methods provide data points automatically. In our 
previous work [27] data points were manually chosen. In this 
paper we used a semi-automatic region growing method for 
providing data points that is a well-known method for seg- 
mentation of medical and non-medical images [28]. To seg- 
ment homogenous regions, semi-automatic region growing 
method is used which requires user interaction to identify 
initial seed point. After determining the seed pixels, the re- 
gion growing algorithm applied using gray scale, spatial in- 
formation and Otsu thresholding method for segmenting the 
region [29]. After applying morphology operation, border of 
the segmented region have been used as initial data points 
for fitting superellipse. In the next step using Levenberg- 
Marquardt algorithm (LMA) [30] for parameter estimation, 
superellipse has been fitted on data points. The algorithm 
diagram can be seen in "Fig. 1,". 

II. Initial Segmentation Using Region Growing Method 

Basically, region growing operates by merging the nearby 
pixels that meet a given homogeneity criterion, starting from 
an initial set of points, known as "seeds". This approach 
offers several advantages over conventional segmentation 
techniques. Instead of identifying boundaries, region growing 
operates always on closed regions in each step of the 
algorithm, and thus avoids further post-processing to recover 
the boundaries of disconnected objects. The algorithm is 
also more stable with respect to noise, and region membership 
can be based on multiple criteria, facilitating the simultaneous 
consideration of several features from the image data, and 
the introduction of eventual a priori knowledge about the 
structures to be segmented. Our proposed semi automatic 
region growing algorithm, needs a seed point that identified 
by operator as shown in "Fig. 2,". Because of almost 
homogeneous left ventricle region in cine MRI, final result of 
region growing algorithm is not dependent to the location of 
seed point and it has minimum effect to the final partial data. 

A. Gray Space map 

The algorithm of region growing is very simple. We 
compute the seed gray level: V, then look for structures which 
have the same gray level than the seed overlapping the seed 
position. At the second iteration, we look for structures 
having a small gray level difference from the seed. In other 
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Figure 1. Our proposed method diagram 
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Figure 2. (a) Manually selection of seed point (b) Visualization of 
Gray Space map 

words we define a set of gray levels from "V-D" to "V+D". (D 
is the difference value). Then we keep those structures which 
overlap the seed position. In each iterations we increase the 
difference D by 1. In this way structures which are closed 
from a spatial AND intensity point of view to the seed are 
highlighted with higher values [32]. In new image if we far 
spatially and from an intensity the point of view from the 
seed, the lower intensity is labeled. Resulted image is Gray 
Space map of the image, "Fig. 2,". An approach is to verify 
that even if the labeled area grows, the homogeneity of the 
area is constant. We try to see statistics of the area during the 
region growing, when the statistics change too much, it means 
that we are introducing heterogeneous areas in our region, so 
we must stop the growing. 

In our case we computed the variation of the Euclidean 
distance between the histograms before and after each 
growing step [29]. We found this variation to be very similar 
to the variation of the region area: if we include large areas in 
our ROI, there is a big area variation, but also, in the same 
time a big statistic variation. "Fig. 3," shows a classical 
histogram after the GS map computation. The ROI is 
highlighted, so it is placed in the higher intensity range. 

B. Segmentation 

First we find the maximum area variation in "Fig. 3," which 
means that from this intensity to is not the Region of Interest 
(ROI). Then we cut the histogram from MAX to intensity 0. 
We have to find the threshold from MAX to the highest 
intensity which separates the uncertainty area from the ROI. 
This is simply done using the well-known Otsu thresholding 
method. This is a parameter free thresholding technique which 
maximizes the inter-class variance. Two types of 
morphological operations are applied to the binary images in 
order to filling this region. We first applied filling algorithm to 
fill the holes. The closing algorithm then applied to have a 
smooth boundary. 

C. Border detection 

The final step for finding initial data points is determining 
of segmented region border. This achieved by applying canny 
edge detection algorithm on segmented region. The results 
are shown in "Fig. 4". 
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Figure 3. Gray Space(GS) map histogram 
III. SUPERELLIPSE FITTING 

Several methods have been proposed for fitting 
superellipses to the images [31]. We used an iterative 
Levenberg-Marquardt algorithm [30] (LMA) for parameter 
estimation. This method provides a numerical solution to the 
problem of minimizing a function, generally nonlinear, over a 
space of parameters of the function. 

A. Distance criteria 

Fitting superellipse on a given set of pixel data using 
least square solution achieved by finding the set of model 
parameters that minimize the sum of the squares of the 
distances between the model curve and given pixel data. The 
notion of distance can be interpreted in various ways and 
Rosin and West [33] investigate several examples. The 
Euclidean distance between a pixel and the point on the 
superellipse closest to it is probably the most suitable 
measure. Unfortunately, finding this distance is generally 
very time consuming and therefore is rarely used. A simpler 
but still effective measure is the algebraic distance given by 



(x - xjcos 0-(y-y c ) 



sin 



(y - yjcos 6-(x-x c ) 



sin 



(1) 



Figure 4. Gray Space(GS) map histogram 

[33]. Tapering parameter, formulated by equation (2), 
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that minimized by another equation Q [ by Levenberg 
Marquardt algorithm. The best fit superellipse is determined 
by finding the parameters which minimize the objective 

functional Q (x, y) 2 . 

One of the main problems with the algebraic distance is 
that it results in different distance estimates for different parts 
of superellipse curves depending on the local curvature of 
the curve. This is because the conventional algebraic distance 
measure treats pixels as individual data points and relations 
between pixels are not exploited. In this paper we propose a 
modification procedure using adaptive data point selection 
and parameter X for more effective convergence of LMA 
algorithm. 

B. Fitting procedure 

Like other numeric minimization algorithms, the 
Levenberg-Marquardt algorithm is an iterative procedure. 
To start a minimization, initial guess must be provided for the 
parameter vector, P [30]. In many cases, an uninformed 
standard guess like p = (l,l,...,l) will work fine; in other cases, 
the algorithm converges only if the initial guess is already 
somewhat close to the final solution. 
So our fitting procedure consists of 3 steps: 

• Determining data points. This step was done by region 
growing method in pervious section. 

• Preparing initial values for Levenberg Marquardt algorithm 

namely X c , y c , a, b, s , Q , t. Our parameter vector, 

therefore should be expressed as p{„ a, b, , , t} 

• Implementation of Levenberg Marquardt algorithm for 
finding superellipse parameters. 

The initial values are given by user. For our propose an 
uninformed standard guesslikeP=(l,l,...,l)will work fine. 
Also data points are given by user. 

After determining data points the Levenberg Marquardt 
algorithm starts fitting superellipse to partial data. First of all 
as a problem definition we have a set of data points of 
independent and dependent variables, (x, y), optimize the 
parameters P of the model curve f(x, P) so that the sum of the 
squares of the deviations 

s(P) = f J [y i -f(x i ,P)] ( 3) 

becomes minimal. 

Levenberg Marquardt algorithm as mentioned is an 
iterative procedure and needs initial values for parameter 
vector p. After preparing theses initial values by user 
algorithm starts its iteration. 

In each iteration step, the parameter vector, P, is replaced 
by a new estimate, P + 5. To determine 5, the 

functions /(x^P + 5) are approximated by their 
linearization. 

f(x i ,p + S)*f(x i ,p + 8) + J i 8. (4) 
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Where 

df(x n p) 

is the gradient (row- vector in this case) of f with respect to p. 
At a minimum of the sum of squares, called S, 
the gradient of S with respect to 5 is 0. Differentiating the 
squares in the definition of S, using the above first-order 

approximation of /(x ; ,P + S), and setting the result leads 
to: 

(J T J)S = J T [y-f(P)]. (6) 

Where / is the Jacobian matrix whose i-th row equals I, and 
where f and y are vectors with ith component 

/(x^P) and y., respectively. This is a set of linear 

equations which can be solved for 5. Levenberg improves 
the algorithm with replacing this equation by a "damped 
version", 

(J T J + lI)S = J T [y-f(P)]. (7) 
Where I is the identity matrix, giving as the increment, 5, 
to the estimated parameter vector, p. The (non-negative) 
damping factor, X, is adjusted at each iterations. If reduction 
of S is rapid, a smaller value can be used, whereas if iteration 
gives insufficient reduction in the residual, X can be increased. 

"Fig. 5," shows the result of applying our proposed 
method on a slice of short axis cardiac cine MRI in whole 
cardiac cycle from end of diastole to end of systole that shows 
robustness of proposed method for segmentation of left 
ventricle in different phases of cardiac cycle. 

"Fig. 6," shows the result of applying our proposed 
method on short axis cardiac MRI images of various patients 
that proves robustness of proposed method for segmentation 
of left ventricle in different shapes. 

IV. Estimating Left Ventricle Volume 

The final step is the calculation of the left ventricle 
volume. Therefore the segmented pixels of all images are 
counted and multiplied by their voxel size. This segmentation 
result, define a binary mask on the input images which 
distinguishes between ventricle and background pixels. In 
the field of radiology the technique of defining the counting 
result of voxels to get a volume is called Simpson Rule. The 
DICOM file format stores the voxel size by means of three 
entries: 

• pixel spacing in x- and y-direction (0.83333 mm x 
0.83333 mm) 

• slice thickness in z-direction (8 mm - 1 .2 mm) 

• slice gap (0 mm - 4mm) 

With this information voxel size is defined as 
pixel _spacing_x * pixel_spacing_y * slice jthickness. The 
volume generated for this approach can be seen in "Fig. 7,". 
"Fig. 8," shows result of implementation of region growing 
algorithm on a not-well-defined left ventricle border short 
axis cardiac MR image that proves undeniable effect of an 
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Figure 6. Result of superellipse fitting on left ventricle in various 
types of cardiac MR images 

expert clinician for providing partial data in these cases. 

V. EXPRIMENTAL RESULTS AND DlSCUTIONS 

Cardiac cine MRI images are used to evaluate left ven- 
tricle volume in this study. Cine MR images provide compre- 
hensible view of a moving organ like heart. The implementa- 
tion of this method is done using MATLAB 7. 12.0. Proposed 
algorithm is tested on 17 data sets of 17 patients. The images 
size was 254x254 pixels. 

To validate the segmentation results and qualitative com- 
parison, we compared obtained results with manual segmen- 
tation performed by a senior radiologist. 

Because of ellipsoid shape of left ventricle, fitting 
superellipse on cardiac cine MRI images has excellent 
accuracy. Our results show less computation and time 
consuming compared to active contour methods for 
segmenting left ventricle. As a quantitative comparison, 
Computational cost reduced in our proposed method, for 
example 1.3 times faster than well-known active contour 
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Figure 8. Region growing algorithm defect. In case of not-well- 
defined LV border, user interaction for providing partial data is 
undeniable 
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method Chan-Vese [34]. One of the limitations which should 
be mentioned is the influence of the papillary muscle, a muscle 
surrounding parts of the left ventricle. The problem is that 
the muscle doesn't fit into the approximated superellipsoidal 
ventricle shape and, although it should be incorporated into 
the ventricle volume, it can't easily be separated from the 
ventricle muscle (myocardium) because both consist of the 
same tissue type and therefore both have a similar grey-value 
on a reconstructed image. In addition radiologists do not 
agree whether the papillary muscle should be included in the 
ventricle or not. 

In this work the papillary muscle will be assumed to be a 
part of the left ventricle and the impact of this assumption on 
the different volume estimation techniques will be 
investigated. 

End-diastolic and end-systolic contours in the images 
were drawn by means of the implemented segmentation tool 
application. Results of the parametric estimation method are 
also required for the evaluation, therefore the end-diastolic 
and end-systolic volume estimation from the patient's 
diagnosis (performed by an operator after the MR scanning 
procedure) were provided. To show that different medical 
operators (and even the same operator at different points of 
time) tend to produce slightly varying estimation results, the 
parametric model was applied a second time on each data set 
by a radiologist. The differing results are of course unintended 
but nevertheless nearly unavoidable. Different levels of 
concentration depending on the number of already processed 
data sets or the tiredness of the operator during the work is 
one possible explanation. Another one is the experience of 
the operator with radiological and (in our case) cardiological 
aspects but also differing interpretations of certain medical 
aspects can lead to varying results. "Fig. 9," illustrates the 
correlation of these two estimations. 

Table I shows all of the provided information of the image 
data sets. Columns"EDV 1 " and"ES V 1 " contain the results of 
the parametric volume estimation while columns"EDV2" 
and"ESV2" contain the older results of the parametric model 
from the actual patient examination. 




Figure 9. Variation in LV segmentation by two clinician experts 

We used four measures to evaluate the segmentation 
results which are: 



2N T 

Similarity index: S t =— p — *100% 

N M +N A 

Where Nj the number of true positive and N M is the 
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Table I. Patients Datasets With Volume Estimation Results From 

SlIPERELLIPSE MODEL FITTING 



P.Num 


EDVl[ml] 


ESVl[ml] 


EDV2[ml] 


ESV2[ml] 


1 


105.9 


37.7 


106.0 


37.9 


2 


173.2 


62.4 


182.0 


67.3 


3 


153.8 


58.5 


153.0 


56.6 


4 


97.1 


62.0 


97.1 


58.1 


5 


172.1 


71.3 


173.4 


66.6 


6 


198.0 


111.8 


202.3 


116.2 


7 


146.3 


101.1 


145.3 


97.5 


8 


125.0 


43.6 


125.0 


43.8 


9 


286.6 


173.3 


288.0 


173.3 


10 


123.4 


56.6 


123.0 


51.2 


11 


83.2 


20.2 


81.3 


20.0 


12 


95.4 


28 


95.2 


32.2 


13 


113.7 


22.0 


114.6 


22.1 


14 


105.5 


18.0 


106.5 


18.2 


15 


98.1 


19.3 


99.1 


20.5 


16 


81.3 


27 


82.1 


29 


17 


77.4 


18.9 


75.2 


16.8 



cardinality of M and N A is the cardinality of A; 



N T 

• Jaccard index: ^' = ~iu 7 m 7 at *100% 

n M + "A + n T p 

• Ratio of correct detection (Sensitivity): 

N T 

77 = — p - * 100% 
P N M 

Np 

Specifity: S E =100-F: Fp = — -^*100% 

Where N p p is the number of false positive. 

Table. II show evaluation of the fitting superellipse on 17 
datasets. Results show that this algorithm can robustly fit a 
superellipse to a MR cardiac image. 



Table II. Evaluation of the fitting results of images for which a manual 

FITTING WAS AVAILABLE 





Volume metric (%) 


SI 


JI 


SE 


SP 


Average 
result 


80.91 


73.28 


81.44 


90.93 



VI. CONCLUTION 

In this work we proposed a semi-automatic method that 
robustly segments and fits superellipse on the left ventricle 
in cardiac cine MRI images. Superellipse can model a large 
variety of natural shapes, including ellipse, rectangles, 
parallelograms, and pinched diamonds, by changing a small 
number of parameters. Several approaches have been 
suggested to determine the parameters of superellipses. We 
used an iterative Levenberg-Marquardt algorithm for this 
approach that needs partial data to fit superellipse on them. 
There are several methods that can provide partial data for 
fitting procedure and a semi-automatic region growing method 
is used for this approach in this work. This method can 
significantly reduce the left ventricle segmentation time that 
is applying manually by radiologists and this is an important 
point because manually segmentation is a long and tedious 
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task, prone to intra- and inter-expert variability, requires about 
20 min per ventricle by an expert clinician. . In this work the 
papillary muscle will be assumed to be a part of the left 
ventricle and the impact of this assumption on the different 
volume estimation techniques will be investigated. We found 
that using this method is less time consuming than active 
contour models and because of ellipsoid shape of left 
ventricle, fitting superellipse on cardiac cine MRI images has 
maximum accuracy. Further researches can work on 
calculation of left ventricle functions like Ejection Fraction or 
stroke Volume from these types of images. They can also 
develop the automatic methods for providing data points 
like ellipse fitting methods. 
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